Journal of Mathematical Biology

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

AO_SCPLOWBSTRACTC_SCPLOWWe study an extension of an environmentally mediated epidemiological model that incorporates direct human-to-human transmission. While the original formulation accounted for environmental exposure, it did not include direct transmission between individuals. Allowing both transmission routes to interact leads to significant qualitative changes in the system dynamics. The analysis reveals multiple dynamical regimes governed by environmental and combined threshold quantities. The stability of the disease-free equilibrium is controlled by an environmental threshold, whereas a combined reproduction number determines the onset of multistability. For certain parameter ranges, endemic equilibria coexist with the disease-free equilibrium, giving rise to backward-type bifurcation behavior and sensitivity to initial conditions. Moreover, the direct transmission rate acts as an organizing parameter by inducing the emergence of an environmental-free equilibrium when exceeding its classical threshold. These results highlight how environmentally coupled transmission mechanisms can generate rich dynamics in low-dimensional models.

We introduce a spectral existence criterion for the evolution of cooperation in the form of the inequality{lambda} maxb > c, where{lambda} max is the leading eigenvalue of an interaction operator encoding population structure, and b and c represent benefit and cost tradeoffs, respectively. Nowaks five rules for the evolution of cooperation correspond to cases in which the cooperation condition reduces to a scalar assortment coefficient. These results follow from the Price equation, which sheds light on a long-standing debate on the role of inclusive fitness and evolutionary dynamics in explaining the evolution of cooperation.

Clavicle fractures often exhibit markedly different clinical outcomes: some patients recover acceptable function despite shortening or displacement, whereas others with apparently similar deformity develop persistent pain, functional loss, or poor healing. To explain this distinction, we propose a minimal nonlinear mechanical model for prognostic analysis of clavicle fractures. The model describes the interaction between fracture-related shortening and compensatory shoulder-girdle posture through a reduced equilibrium equation incorporating stiffness, geometric nonlinearity, and shortening-posture coupling. Within this framework, we analyze equilibrium branches, local stability, and the emergence of critical thresholds. We show that post-fracture destabilization can be interpreted as a fold bifurcation, while more complex parameter dependence gives rise to cusp-type structures and multistability. These bifurcation mechanisms provide a mathematical explanation for sudden deterioration after injury or treatment, as well as for strong inter-individual variability. We further introduce an optimization principle based on a utility functional to guide treatment planning. The analysis predicts that the optimal safe correction should lie strictly below the bifurcation threshold, thereby generating a natural safety margin. Although the model is simplified and has not yet been calibrated against patient data, it nevertheless provides a theoretical framework for understanding why fracture prognosis may deteriorate abruptly near critical mechanical conditions and offers a dynamical-systems interpretation of empirical treatment thresholds used in clinical practice.

The analysis of large gene and metabolic networks is often hindered by unknown biochemical parameters and the nonlinear nature of classical S-system models. To address this, we introduce a framework based on combinatorial toric geometry computed with tools such as Normaliz, SageMath, it is worth mentioning this technique in not restrictive to integer vectors, there exists a natural extension to real geometries. Unlike traditional approaches, which rely on parameter dependent fixed points, our method constructs a Topological Environment derived from the dual space of kinetic orders, leading to what we call orthogonal enzyme kinetics. Within this topological setting, fixed points are computed on the algebraic torus, enabling the transformation of nonlinear dynamics into linear forms. Importantly, these fixed points are independent of kinetic parameters and depend only on network topology and interaction signs. Applying this methodology to gene circuits involved in circadian rhythms, we reproduce previously reported oscillatory physiologies.

Recent analyses of the effects that organismal genealogies or pedigrees of populations have on times to common ancestry for samples of genetic data are extended to cases of population subdivision and migration. Traditional coalescent models marginalize over pedigrees. A finding of a pedigree effect implies that data analysis and interpretation should not be based on the corresponding traditional coalescent model but rather on a coalescent model obtained by conditioning on the pedigree. We apply a straightforward test based on the distribution of pairwise coalescence times to four previously described scenarios of subdivision and migration. These scenarios are defined by the relative magnitudes of four parameters: the number of the local populations or demes, the deme size, the migration fraction, and the probability that migration can occur at all. We find pedigree effects in three scenarios. In two, the effect is weak if the deme size is large. The one scenario without a pedigree effect corresponds to the well known structured-coalescent model. The one scenario with a persistent pedigree effect even in the limit as the deme size tends to infinity involves long periods without gene flow interrupted by pulses of migration. We illustrate our results using simulations and numerical analysis.

Biological tissues grow at rates that depend on the geometry of the supporting tissue substrate. In this study, we present a novel discrete mathematical model for simulating biological tissue growth in a range of geometries. The discrete model is deterministic and tracks the evolution of the tissue interface by representing it as a chain of individual cells that interact mechanically and simultaneously generate new tissue material. To describe the collective behaviour of cells, we derive a continuum limit description of the discrete model leading to a reaction-diffusion partial differential equation governing the evolution of cell density along the evolving interface. In the continuum limit, the mechanical properties of discrete cells are directly linked to their collective diffusivity, and spatial constraints introduce curvature dependence that is not explicitly incorporated in the discrete model. Numerical simulations of both the discrete and continuum models reproduce the smoothing behaviour observed experimentally with minimal discrepancies between the models. The discrete model offers further individual-level details, including cell trajectory data, for any restoring force law and initial geometry. Where applicable, we discuss how the discrete model and its continuum description can be used to interpret existing experimental observations.

Biological organisation is inherently multi-level: molecular processes, membrane dynamics, cellular geometry and tissue context reciprocally constrain one another, often through boundary-mediated feedback. A recurring theme in theoretical biology is that such organisation is not well captured by models that assume a fixed repertoire of variables and a pre-given state space: what counts as a relevant state description can depend on organisational context and history. The principle of biological relativity further sharpens the same challenge from a different angle, emphasising that no level is causally privileged and that cross-level feedback can close into circular causality. These lines of work motivates for a structural multi-level semantics for modeling the biological pathways. We introduce a constraint-based semantic framework that distinguishes an evolving organisational scaffold--the admissible multi-level patterns and interfaces--from the pathways that traverse and coordinate them. This separation yields mathematical, loop-level diagnostics for boundary-driven circular causality: it identifies when organisational trajectories induce persistent reparameterisations of local state descriptions, and it classifies cyclic regimes into reversible loops, stable history-dependent loops, and unique (rare) organisational reconfigurations. The framework is accompanied by a systematic crosswalk to mainstream causal, dynamical and computational approaches, clarifying what is gained when interfaces and local-global consistency are treated as semantic, rather than purely parametric, structure. We demonstrate the approach on a canonical excitable-cell exemplar by modelling a single Hodgkin spike as a cross-level interface loop coupling membrane, molecular and cellular constraints. Without re-deriving Hodgkin-Huxley kinetics, the resulting diagnostics provide an explicit semantics for boundary-mediated feedback and spike-induced history dependence, including when cyclic activity imprints persistent changes in effective excitability. Together, the case study and comparisons position constraint semantics as a practical mathematical layer for multi-level biological organisation: compatible with existing mechanistic models, yet designed to expose circular causal closure and organisation-dependent state descriptions that standard formalisms typically leave implicit. AMS subject classifications92C30, 92C46, 92B05, 55U10, 55R10

Vertical transmission of an infectious disease from parent to offspring is a common transmission route in many systems. Here we investigate the dynamics of a pathogen with both horizontal and vertical transmission within a spatially structured population. We introduce a lattice model with a pair approximation that includes both local and global transmission and reproduction. We find that vertical transmission can determine pathogen invasion and reduce the horizontal transmission rate required for invasion. When the majority of transmission and reproduction is local, vertical transmission can destabilise a host population to cause limit cycles. Given the advantages of a pathogen having both horizontal and vertical transmission routes, we extend the model to investigate the likelihood a mutant strain with both transmission modes will outcompete a resident strain with only horizontal transmission. When there is no trade-off the mutant always invades and when there is a trade-off with horizontal transmission, the mutant emerges when the cost to the horizontal transmission rate is not too large. Depending on how the mutant appears within the host population, it may have an initial advantage over the resident strain even if it cannot outcompete in the long-term. Our work demonstrates the potential importance of vertical transmission within host-pathogen dynamics.

Cell-cell adhesion is a key regulator of cancer invasion. In this work, we extend a pre-existing individual based cancer invasion model by introducing a stochastic representation of N-cadherin-mediated adhesion, where the lifetime of a cell-cell bond depends on the pulling force acting on the bond. Using experimental data, we derive expressions for the mean and standard deviation of N-cadherin bond lifetimes and fit them to Gamma distributions, enabling their treatment as force-dependent random variables. These distributions are then used to modify the diffusion coefficient of mesenchymal cancer cells. The model predicts reduced random motility with increasing adhesion and incorporates a dynamic transition between catch- and slip-bond behaviour. Along with this model for cell motility, we propose a preliminary physical framework, that can be used to model pattern formation as a result of the new adhesion mechanic.

Imagine you walk in a plane. You move by making a step of a certain length per time interval in a chosen direction. Repeating this process by randomly sampling step length and turning angle defines a two-dimensional random walk in what we call comoving frame coordinates. This is precisely how Ross and Pearson proposed to model the movements of organisms more than a century ago. Decades later their concept was generalised by including persistence leading to a correlated random walk, which became a popular model in Movement Ecology. In contrast, Langevin equations describing cell migration and used in active matter theory are typically formulated by position and velocity in a fixed Cartesian frame. In this article, we explore the transformation of stochastic Langevin dynamics from the Cartesian into the comoving frame. We show that the Ornstein-Uhlenbeck process for the Cartesian velocity of a walker can be transformed exactly into a stochastic process that is defined self-consistently in the comoving frame, thereby profoundly generalising correlated random walk models. This approach yields a general conceptual framework how to transform stochastic processes from the Cartesian into the comoving frame. Our theory paves the way to derive, invent and explore novel stochastic processes in the comoving frame for modelling the movements of organisms. It can also be applied to design novel stochastic dynamics for autonomously moving robots and drones.

The immune system was historically defined as a system that provides protection from pathogens. Numerous models have been developed to understand how immunity faces a complex world of microbes that includes pathogens and symbionts, as well as cells of our own self that may develop tumors. Based on the classical assumption that survival depends on internal homeostasis, we have developed a formal model of homeostasis for a host interacting with microbes and self. We propose that such a model must include two fundamental functions: a function that counters change (including tissue repair), and a function that counters the agent of change (such as "immunity" to microbes or self). We show that this elementary model is sufficient to generate symbiosis, and that symbiosis is an emergent property of the host-microbe relationship that does not require the microbe or the host to express "traits of symbiosis". We suggest that the conditions leading to symbiosis contribute to eukaryotic evolution and ontogeny. This model may be further applied to symbiotic interactions between organisms and non-microbial or non-cellular agents of change.

We present an approach to analysing cell homeostasis using a bond graph modelling approach that ensures that the conservation laws of physics (conservation of mass, charge, and energy, respectively) are satisfied for the interdependent biochemical, electrical, mechanical, and thermal energy storage mechanisms operating within the cell. We apply the bond graph approach to several cell membrane transport mechanisms and then consider how physics constrains intracellular electrolyte homeostasis for enterocytes (the epithelial absorptive cells of the gut). The model includes the electrogenic sodium-potassium ATPase pump (NKA), the glucose transporter (GLUT2), and an inwardly rectifying potassium channel, all in the basolateral membrane, and the electrogenic sodium-driven glucose transporter (SGLT1) in the apical membrane. Glycolysis converts the imported glucose to ATP to drive NKA. For specified levels of sodium, potassium, and glucose in the blood, the model demonstrates how enterocytes absorb sodium and glucose from the gut and transfer glucose to the blood while maintaining the membrane potential and homeostasis of intracellular sodium and potassium. The Gibbs free energy available from the ATP hydrolysis ensures that the cell operates as a sodium battery with a high external to internal ratio of sodium concentration in order to provide the energy for many other cellular transport processes. We show that the 3:2 stoichiometry of Na+/K+ exchange in NKA, coupled with 2:1 Na+/glucose cotransport in SGLT1, a 1:2:2 ratio between glucose consumption and ATP and water production in glycolysis, and K+ and glucose efflux through Kir and GLUT2, respectively, provides a balanced system that maintains homeostasis of intracellular Na+, K+, glucose, ATP and water, and homeostasis of the membrane potential, under varying levels of transport of glucose from the gut to the blood. We also show how the flux expressions for SLC transporters, ATPase pumps and ion channels can all be expressed in a consistent and thermodynamically valid way.

The basic reproduction number of an infectious disease is known to depend on the structure of contacts between individuals in a population. This relationship has been explored mathematically through two well-known models: one which depends on a matrix of contact rates between different demographic groups, and another which depends on the variability of contact rates over the population. Here we introduce a model that combines and generalises these two approaches. We derive a formula for the basic reproduction number and validate it through comparisons to simulated outbreaks. Applying this method to contact survey data collected in Belgium between 2020 and 2022, we find that our model produces higher estimates of the basic reproduction number and larger relative changes over periods when social contact behaviour was changing during the COVID-19 pandemic. Our analysis suggests some practical considerations when using contact data in models of infectious disease transmission.

A common assumption in molecular evolution is the fixed selection nature of a mutation, for instance, a neutral mutation is selectively neutral for all individuals who carry the mutation, and so forth a deleterious or beneficial mutation. Our recent work challenged this presumption, postulating that individuals with a specific mutation exhibit a fluctuation in fitness, short for FSI (fluctuating selection among individuals). Moreover, an intriguing phenomenon called selection-duality emerges, that is, a slightly beneficial mutation could be a negative selection (the substitution rate less than the mutation rate). It appears that selection-duality is bounded: the low-bound is the generic neutrality where the mutation is neutral by the means of fitness on average, while the up-bound is the substitution neutrality where the substitution rate equals to the mutation rate. In this paper, we conducted a thorough theoretical analysis to evaluate how many generations needed for a selection-duality mutation to be fixed in a finite population. A striking finding is that the mean fixation time of a selection-duality mutant, including the generic neutrality and the substitution neutrality, is approximately identical, which is considerably shorter than the case of strict neutrality without FSI. One may further envisage that the fast-fixation nature of selection-duality mutations could result in a considerable genetic reduction at linked sites.

We introduce a mechanistic, nonlocal tumour-growth model designed specifically to capture explosive dynamics that are not adequately explained by standard logistic reaction-diffusion descriptions. The motivation is empirical: the universal scaling law reported in [1] provides compelling cross-sectional evidence of superlinear tumour activity versus tumour burden, but as a phenomenological relationship it does not by itself supply a dynamical mechanism, nor does it rigorously describe how explosive growth emerges, how fast it develops, or how spatial interactions and tissue boundaries influence it. Our model addresses this gap by incorporating nonlocal proliferative feedback--cells respond to a spatially aggregated neighbourhood signal--and a singular, Kawarada-type acceleration that produces "quenching": tumour density stays bounded while the proliferative drive becomes unbounded as the aggregated signal approaches a critical threshold. This offers a concrete mechanistic route to explosive escalation consistent with physical boundedness. We analyse the model under no-flux (Neumann) boundary conditions, appropriate for reflecting tissue interfaces. In the spatially homogeneous setting we prove finite-time onset of the explosive regime and obtain explicit rates for how rapidly it is approached. For spatially heterogeneous perturbations we derive a transparent spectral stability theory showing how the interaction kernel selects spatial scales and how the singular acceleration tightens stability margins as the explosive threshold is approached. These results provide interpretable links between nonlocal interaction structure, boundary effects, and the emergence of rapid growth. Finally, to connect mechanism to data in the spirit of [1], we embed the model in a Bayesian inference framework that treats the interaction kernel and the acceleration strength as unknown and learned from tumour-growth observations. This enables uncertainty-aware estimation of explosive onset times, escalation rates, and stability margins, while positioning the scaling law of [1] as an observable signature that our mechanistic model can explain and quantify rather than merely fit.

In this research, we create a new fractional-order SEIHRD framework to examine how the Nipah virus moves from one species to another (zoonotic spillover) and how it later spreads throughout a community (via contact with one another) or in a hospital or isolation situation (via entering into a hospital or being placed under quarantine). We used the fractional-derivative formulation of the SEIHRD model to demonstrate memory-based effects related to the progression of an infection and also reflect time-distributed effects associated with surveillance and control measures placed on an infected patient. We first demonstrated that the basic epidemiologic properties of the model were consistent by showing that the solutions of the SEIHRD differential equations will always yield positive and bounded solutions within biologically relevant parameter ranges. We then established the well-posedness of this model by transforming the SEIHRD differential equations into an equivalent integral operator and applying various fixed-point arguments to demonstrate that there will always be unique solution(s) to the SEIHRD differential equations. To evaluate the threshold parameter for the transmission of Nipah virus within a given population we calculated the threshold level through the next generation method to determine the expected number of secondary infections from a new or chronically infected host. One of the main contributions of this work is to include an analysis of the robustness of a given solution to all potential perturbations (i.e., Ulam-Hyers and generalized Ulam-Hyers stability). In addition, we provide analytic results guaranteeing that small perturbations due to approximate modeling, numerical approximation (discretization), or the lack of data fidelity will produce controlled deviations in the solutions. To finish this project, we perform a global sensitivity analysis on uncertain coefficients to evaluate their contribution to the uncertainty of each coefficient and to find out the coefficients that most strongly influence major outcome metrics. This will allow us to develop a priority order for prioritizing spillover control (reduction of human contact and/or isolation), contact reduction, and expenditure of resources towards isolation-related interventions. The resulting framework converts fractional epidemic modeling from a descriptive simulation to a replicable method with robustly defined behavior and equal response prediction.

Cancer-associated fibroblasts (CAFs) are a central component of the tumor microenvironment that facilitate a supportive niche for cancer progression and metastasis. Experimental evidence suggests that CAFs can facilitate estrogen-independent tumor growth, thereby reducing the efficacy of anti-hormonal therapies. Understanding and quantifying the complex interactions between tumor cells, hormonal signalling, and the microenvironment are crucial for designing more effective and individualized treatment strategies. We propose a mathematical framework to explore the influence of CAFs on ER+ breast cancer progression and to evaluate strategies to mitigate their impact. We develop a system of nonlinear ordinary differential equations that substantiates the experimental observations by providing a mechanistic basis for the role of CAFs in regulating estrogen-independent growth dynamics. We then employ optimal control theory to evaluate distinct therapeutic approaches involving monotherapy or combinations of: (i) inhibition of tumor-to-CAF signaling, (ii) inhibition of CAF-to-tumor proliferative signaling, and (iii) endocrine therapy. Taken together, our results demonstrate that CAF-targeted strategies can enhance treatment efficacy across various estrogen dosing regimens. Our study provides new insights into the potential of CAF as a therapeutic target that could help to improve existing approaches for endocrine therapies.

Organisms must manage a trade-off between robustness and flexibility as they enact adaptive behaviors. One way organisms achieve this is by navigating a network of quasi-stable behavioral states. Evidence for such behavioral states has been observed in many organisms, and new methods for detecting these states have taken on a prominent research focus. Although dynamical models demonstrating behavioral switching have been developed significantly over the past few decades, theories of the similarities and differences among these models, necessary for advancing empirical modeling, have not yet been fully elaborated. Here, we consider behavioral switching in two different classes of dynamical models of the forward-reversal behavioral transition in C. elegans. We first show how fundamentally different models can give rise to similar phenomena under noisy conditions. We then analyze the deterministic aspects of these models to expand on their differences, clarifying the theoretical relationship between them. Finally, we demonstrate how sequence models can be further extended to incorporate dwell times for behavioral states. Our work contributes toward a broader theoretical understanding of behavioral switching in adaptive systems.

Cytokine-mediated communication is a central mechanism by which immune cells coordinate activation, differentiation and proliferation. While mechanistic reaction-diffusion models provide detailed descriptions of cytokine secretion and uptake at the cellular scale, their computational cost limits their applicability to large and densely packed cell populations. Previously employed approximations of cytokine diffusion fields rely on assumptions that neglect the influence of cellular geometry and volume exclusion. In this work, we study a macroscopic description of cytokine diffusion and reaction dynamics based on homogenization techniques, rigorously linking microscopic reaction-diffusion formulations to effective continuum models. The resulting homogenized equations replace discrete responder cells with a continuous density, while retaining essential features of cellular uptake and excluded-volume effects. Further, we show that in regimes with approximate radial symmetry, classical Yukawa-type solutions emerge as limiting cases of the homogenized model, provided appropriate correction factors are included. Overall, our approach allows efficient multiscale modeling of cytokine signaling in complex immune-cell environments.